FINAL REPORT 


for 

Computer Simulation of Surface and Film Processes 

November 1 » 1 984 


Cooperative Agreement No.: NCC 2-125 

Period of Award: May 1, 1981 - April 30, 1984 

Principal Investigator: Professor W. A. Tiller 

Senior Investigator: Dr. M. T. Halicioglu 

(NASA-CR-176966) ; COMPUTER SIMULATION OF 

SURFACE. AND FILE ■ PROCESSES Final Report, 1 : 

May 1981 - 30=Apr. 1984 (Stanford Oniv.) 

■jc' D CSCL 20L 

V G3/76 


N 8 6-3 1377 

Unclas 

43016 


SU-DMS-85-R-1 


Stanford/NASA Ames Joint Institute for Surface and 
Microstructure Research 

Department of Materials Science and Engineering 
Stanford University, Stanford, CA 94305 



CHAPTER I 


Introduction 

A detailed knowledge of the structure at an atomistic level is very 
important for understanding of various processes which take place at 
surfaces and interfaces. The exact geometry of different atomic configura- 
tions is needed for a correct interpretation of various experimental obser- 
vations. In spite of the sophisticated experimental techniques available 
for examining surface structure and surface chemistry, many structure- 
related characteristics (e.g., surface reconstruction, catalysis) remain 
unresolved. Specific surface processes like catalysis, diffusion, roughen- 
ing, wetting, adhesion, nucleation, crystal growth, crack propagation, 
corrosion, etc. , are so complex that detailed understanding of surface 
structure is essential for understanding the key elements involved in the 
specific process. Unfortunately, conventional theoretical techniques have 
great difficulty in coping with the many "nested" phenomena involved in 
these processes at a detailed level. 

Computer simulation techniques, on the other hand, employing advanced 
numerical methods and well-known statistical approaches, are now beginning 
to contribute significantly to the microscopic understanding of surfaces 
and ad-layer technology as well as to a better understanding of the atomic 
mechanisms involved in various materials surface-related processes [1,2]. 

The main issue in computer simulation is, of course, the ability to 
generate a mathematical description or a "model" that will accurately 
reflect what is happening in the real world. Prior to the advent of 
‘ computer technology, applied mathematics provided a few situations which 
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could be described by analytically solvable equations. Unfortunately, 
important aspects of surface-related processes and materials behavior 
usually are complex phenomena which tend to be resistant to analytical 
treatment and even to definitive studies by laboratory experimental means. 
During the past decades, complex processes have been studied, with increas- 
ing frequency, using "computer simulation experiments." The results have 
been particularly useful and illuminating [1-3]. These simulation studies 
have motivated not only improved laboratory experiment designs but also 
have made possible improved data interpretation techniques. 

The most important objective of this project was to acquire an atomic 
level information for different surface-related phenomena. This was 
accomplished using basically three different simulation techniques [4]. 

All the investigations which were carried out employed in one way or 
another a computer simulation technique based on atomistic level considera- 
tions. In general, three types of simulation methods are being used for 
modeling systems with discrete particles that interact via well defined 
potential functions. 

(i) Molecular dynamics: this is a general method for solving the 

classical equations of motion of a model system. It provides 
time evolution of a system of many particles. One can obtain 
pertinent phase space trajectories, therefore, any physical 
property can, in principle, be calculated if the system 
Hamiltonian is known. Given the initial position and velocity 
vectors for each particle, the dynamical history of the 
assembly is generated in the computer. The simultaneous equa- 
tions of motion are solved by any of several techniques which 
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involve extrapolating the atomic trajectories over successive 
short time intervals. After each extrapolation the forces on 
all the particles are recalculated at their new positions, and 
the resulting accelerations used in the next extrapolation. 

(ii) Monte-Carlo: in general, the term Monte-Carlo refers to the 

Markov chain ensemble averaging technique introduced by 
Metropolis et al. It is purely a stochastic method to model 
equilibrium properties of a system. Accordingly, Monte-Carlo 
methods cannot simulate atomic motions in real time, but gener- 
ate a large number of configurations (i.e., snapshots of N- 
particle positions representative of an ensemble at T). 
Equation of state data are obtained by appropriate averages 
over these generated configurations. 

(iii) Molecular statics: this method, in principle, can only provide 

properties of a system at T = 0°K. Its widespread use stems 
from its small demand of computer time and its ability to 
handle large systems. Molecular statics can be regarded as a 
minimization procedure. In general, the objective is to find 
the energetically most favorable configuration. It may be 
thought of as molecular dynamics with a damping force which 
progressively draws energy out of the system until it arrives 
at a configuration of stable static equilibrium (i.e., a mini- 
mum of the potential energy function). Frequently, molecular 
statics results are employed as starting configurations (i.e., 
as an input) in the Monte-Carlo and molecular dynamics calcula- 
tions. 
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All three of these techniques are based on a well-defined potential 
function describing interactions among particles in the system. 

In general, for a system of N particles, the total potential energy 
may be expanded as [ 5 ] : 



where, u(r,,r.), u(r.,r.,r t ), ..., u(r.,r., ..., r ) are 
i j i j k x j n 

two, three, and n-body potentials, respectively. The position of the i^ 
particle is denoted by r^ 

Clearly, the most important term in this expansion is the first term 
involving two-body interactions. Therefore, in the majority of the 
atomistic calculations made to date, only pairwise additive potentials have 
been used. This provides great simplification in the analytical formalism 
as well as in numerical computations. 

In this project, depending on the type of the system under considera- 
tion, the total potential energy $ was calculated either considering 
"only" a two-body term, or for more quantitative results it was calculated 
as a sum of two- and three-body interactions neglecting four- and higher- 
body terms in eq. (1). Because this expression (eq. (1)) had to be used in 

lengthy machine computations, u(r. ,r.) and u(r. ,r.,r.) were 

i j 1 J K 

' - chosen with the simplest possible functional forms. In this study, 
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therefore, the two-body part was represented by a.Mie-type potential which 
is given by: 


ij ij 


( 2 ) 


" lth r ij - l r i - r ji 


’ r 0 


the energy at r^ = r Q . 
sive and attractive terms 
hand, was expressed as: 


denotes the equilibrium separation and e is 
The exponents m and n account for the repul- 
respectively. The three-body part, on the other 
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(3) 


where, the summation includes all triple multipole interactions resulting 

from the expansion of the third-order interaction energy for three atoms. 

Each term in the summation is expressed as the product of a geometrical 

factor G(r.,r.,r,) which depends on the relative positions of the 
1 J 

three atomic nuclei and an interaction constant which depends only on the 
atomic species involved in the interaction. The functional forms of 
GCr^r^r^) for several multiple interactions have been obtained by 
Doren and Zucker [6]. Here, we consider only the triple-dipole interaction, 
which has been shown to be the dominant contribution [7]. This term was 
first obtained for closed shell atoms as: 
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( 5 ) 
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where 0^, 0_. , 0^ and r„ , represent the angle and the sides of 

the triangle formed by the three particles i, j and k. 

In general, model systems interacting only via two-body potentials 
have been used for parametrical analyses, or to obtain semi-qualitative 
results for resolving simple mechanisms at atomistic levels. While these 
two-body model potentials are able to reproduce characteristics of some 
systems (i.e. , solid rare gases, fee crystals, etc.), they are unable to 
simulate other systems with somewhat more complex structures (e.g., 
diamond, graphite structures). Therefore, for systems with more involved 
structures three-body interactions are included by combining eqs. (1) 
through (5). Accordingly, the total potential energy was expressed as: 


§ 


■ 1 1 ? <^> 

1 J ij ij 


, Z(1 + 3 Cos 0, Cos 0 O Cos 0 O 

•jr I I I — 5 S 1 

■ 1 j k Cry • r lk • r k > 


( 6 ) 


This equation has been used earlier for other analyses, when the 
effect of three-body interactions on the structural characteristics of 
small clusters [8] and on various crystalline solids [5] were investi- 
gated. In these analyses the importance of the many-body effect has been 
clearly demonstrated. The function $ containing three-body interactions 
is able to provide stability regions for many different types of 
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crystalline materials [5]. In this project, to investigate the multi-body 
effect further, in Chapter II we analyze the effect of three-body forces on 
the vibrational frequencies of triatomic clusters. In another study, which 
is presented in Chapter III, the multilayer relaxation phenomena for low 
index planes of an fee crystal was analyzed also as a function of the 
three-body interactions. 

For the simulation of real systems, the energy parameters of eq. (6) 
are needed. Evaluation of these parameters is, in general, a cumbersome 
procedure which is outlined in Chapter IV along with the numerical values 
of the parameters for some selected compounds. In Chapter V various 
surface properties for Si and SiC systems are calculated. 

The rest of the studies in this report are related to materials 
applications that involve analyses of responses of materials to external 
forces. In Chapter VI, results obtained from static simulation calcula- 
tions for slip formation are presented, while in Chapter VII more elaborate 
molecular dynamics calculations on the propagation of cracks in two- 
dimensional systems are outlined. 
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CHAPTER II 


The Effect of Three-body Forces on the Vibrational 
Frequencies of Triatomic Clusters 

In general, experimental information about the energetics and 
structural characteristics of very small "isolated" clusters are obtained 
from various spectroscopic measurements [9-11]. Formal relationships 
between the observed spectral lines and the interatomic forces operational 
among the particles in the cluster are well established. From a knowledge 
of. the interatomic potentials, for example, vibrational frequencies associ- 
ated with a given cluster can be calculated [12-15], For this purpose, we 
employed a normal mode procedure. After an energetically stable configura- 
tion for the cluster is found, the total potential energy is expanded about 
this minimum from which the force constants are determined [14,15]. A 
diagonalization of the force constant matrix produces the desired eigen 
frequencies. This procedure is for the T = 0°K case and is based on the 
harmonic approximation. 

For a system of three particles, eq. (6) may be put in a dimensionless 
form as: 


<£* = 


3 3 . 

I I l(^) 

i J ij 
i>j 


12 


2(75-) 6 } + z* 

ij 


1+3 Cos 


Cos 02 Cos 0^ 


(7) 
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r* )' 
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where the reduced quantities are defined by 
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Z r.. 

Z* = — . and r* . = -2J- . 

9 1J r n 

£r 0 0 

Here, we also set m = 12 and n = 6 which reduces the two-body part of 
the potential to a Lennard-Jones function. This assumption simplifies the 
comparison of the present results with the earlier investigations [8,16-19] 
where the energetics and structural stabilities of clusters have been 
analyzed. The reduced potential energy, £>*, when varied as a function of 
Z*, produces either a linear or an equilateral triangular shape to be the 
energetically more stable form [8]. This is demonstrated in Fig. 1 for 
three different Z*-values. In Fig. 1, the curves represent minimized 
total energies versus the angle 0 which was varied from 50° to 180° so as 
to include both configurations, equilateral triangle and linear. For 
Z* = 0.2, no minimum is discernible for the linear case (9 = 180°). For 
Z* = 0.6 and 1.0, there are noticeable minima at 0 = 180° which are more 
shallow than the equilateral triangular configuration (9 = 60°). Other 
calculations [18-22] support this trend indicating that the energetically 
favorable structures would be either linear or "nearly" equilateral 
triangular. 

Fundamental frequencies for the equilateral triangular and the linear 
configurations were estimated from corresponding characteristic equations 
based on the normal mode procedure. Each configuration exhibits a spectrum 
with frequencies and Figures 2a and 2b display calculated 

"reduced" frequencies w* as functions of the three-body intensity param- 
eter Z*, for the triangular and linear cases, respectively. Reduced 
frequencies were calculated as oj* = cj^/gOq where u)q denotes the 
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corresponding diatomic normal mode frequency. For homonuclear cases, based 
on eq. (2), it is given by 

io 0 = (2mns/M) ^ 2 /r 0 (8) 

where M represents the atomic mass. For the linear configuration, the 
three main peaks corresponding to bending, symmetric and asymmetric vibra- 
tions were well separated. The bending mode exhibits the lowest frequency 
which is expected due to the shallow minimum at 6 = 180° shown in Fig. 1. 
The peak corresponding to the asymmetric vibration, on the other hand, was 
found to be in the high frequency region. For the equilateral triangular 
case the asymmetric and bending modes become degenerate, as anticipated, 
and are located in a lower frequency region than the peak for the 
symmetrical vibration. These trends, for either configuration, remain 
unaltered for the range of Z* (from 0.0 to 1.0) considered in this 
investigation. Calculated peak positions for the equilateral triangular 
configuration are consistent with results reported by Etters et al. [18] 
for the rare gas clusters. Their results represent a special case in our 
approach with Z* = 0. 

The effect of Z* exerted on the peak positions for the linear and 
triangular configurations was found to be in opposite directions. For 
increasing three-body intensities, the two main vibrational peaks for the 
equilateral triangular shape shift to the lower frequency region, while the 
three peaks of the linear configuration shift to the higher frequency 
domain. The largest shift was found to be exhibited by the symmetric 
vibrational mode for the equilateral triangular case. Probably, the most 
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interesting region of Z* values lies between 0.66 and 0.70 where the two 
configurations may coexist in appreciable concentrations. For Z* = 0.68 
the equilateral triangular and linear shapes (with m = 12 and n = 6) 
become energetically degenerate. Accordingly, for systems with homonuclear 
triatomic clusters, one would anticipate a mixture of these species at 
roughly equal concentrations. Comparison of the calculated results with 
experiments for this special situation is difficult, not only because of 
the parametrical nature of the present study, but also because of various 
complications involved in the interpretatin of the spectroscopic measure- 
ments. Recently, Richtmeier et al. [21] calculated vibrational frequencies 
for group IB trimers. For triangular configurations (not necessarily equi- 
lateral) they found that the symmetric stretch mode exhibited the highest 
frequency, while the bending and asymmetric vibrational frequencies were 
lower and relatively close to each other. Basically, these results are 
consistent with the present calculations based on the chosen parameters. 
Some effects due to different m and n values on the relative position 
of vibrational peaks are expected but they were not analyzed in this 
investigation. 

Conclusion 

The calculated normal mode frequencies of homonuclear triatomic 
clusters were found to be affected considerably by the intensity of the 
three-body forces operational among the particles. For increasing intens- 
ity of the three-body forces, the two vibrational peaks for the equilateral 
triangular configuration shift to lower frequencies while the three peaks 
- for the linear shape shift to the higher frequency domain. For homonuclear 
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triatomic clusters with Z* = 0,68, the two species, i.e., the equilateral 
triangular and linear configurations, may coexist. 
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Figure Captions 


Fig. 1 Variation of the reduced potential energy, as a function of 9 
for three different Z*-values. 

Fig. 2 Reduced frequencies versus Z*-values; (a) for the equilateral 
triangular configuration; (b) for the linear configuration. 
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CHAPTER III 


Multilayer Relaxation Calculations 
for Low Index Planes of an FCC Crystal 

In many studies related to surfaces the relaxation and the reconstruc- 
tion of the exposed region are very important. Numerous experimental and 
theoretical works have been generated to explore and understand the atomic 
nature of the surface reconstruction phenomenon. Recent experimental find- 
ings indicate that in many metal surfaces a multilayer relaxation takes 
place. Surface region interlayer spacings differ from the bulk value 
considerably. In general, the most pronounced relaxation occurs in the 
spacing between the first and second atomic layers, d^ 5 which is found to 
be less than the corresponding bulk value [23-31], These experimental 
results have been predicted first by Finnis and Heine [32] using a single 
layer relaxation approach. More recently Landman, Hill and Mostoller [33] 
calculated the reconstruction in metal surfaces using a multilayer relaxa- 
tion procedure and obtained results consistent with experiments. In both 
of these approaches, calculations were performed employing surface elec- 
tronic charge densities based on quantum mechanical considerations. Other 
theoretical calculations using only semi-empirical pair potentials, how- 
ever, predict d^ to be larger than the bulk value which is in contradic- 
tion with experimental results. This shortcoming of the semi-empirical 
pair interaction models which are being used in many computer simulation 
calculations (because of their functional simplicity) inhibits their usage 
in various surface region modeling studies. 
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In this investigation, in order to improve the applicability of the 
semi -empirical potential models in calculations related to the surface 
structures, a potential energy function comprising two-body and three-body 
interactions was taken into consideration. As model systems (100), (110) 
and (111) index planes of an fee crystal were employed. Three-body forces 
were found to be extremely important in the multilayer relaxation of these 
surfaces. The results of calculations showed a trend which was in good 
agreement with experimental findings. 

Throughout this investigation the surface energy a per atom was 
calculated from: 

M 

a = l (e. - e Q ) (9) 

Z 

where e^ and e^ denote the total potential energy for an atom located 
in the X'th surface layer (from the top) and of an atom located in the bulk 
of the system, respectively. M is the total number of the surface layers 
considered in the calculation. 

The cut-off radius, R , was taken to be approximately 5rQ, in all 
cases. In the calculations of eg (the energy per atom in the bulk) the 
stability condition for the crystal was taken into account by assuming: 

5e 0 

W 2 - 0 <10 > 

with V denoting the total volume. For different Z, the values of e^ 
were calculated considering eq. (6) based on the condition imposed by 
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eq. (10). This approach has been employed in reference [5] in calculating 
stabilities of ideal crystals. 

Results and Discussions: 

In numerical calculations, the units for the energy and distances were 
reduced by e and Tq (of eq. (2)), respectively. For the multilayer 
relaxation procedure the surface and bulk energies (i.e., e^ and Cq) were 
calculated for different Z values ranging from 0.0 to 0.4. Three dif- 
ferent surface planes (100), (110) and (111) for an fee crystal were in- 
cluded in the calculations. Relaxations were performed for the top-most 
five surface layers by minimizing the surface energy a of eq. (9) with 
respect to the positions of the atomic layers. In this minimization pro- 
cedure, for every Z value, the bulk energy e^ was calculated consider- 
ing eqs. (6) and (10), simultaneously. The minimization involves a layer- 
by-layer relaxation which was carried out by varying the positions of the' 
layers only in the perpendicular direction (z-direction). However, atomic 
arrangements within each layer (i.e., x-y planes) were left unaltered. 

In Fig. 1, the energies calculated for three different surface 
structures and for the bulk were plotted versus the Z value. For all 
three cases, the surface energies were found to be monotonically decreasing 
functions of Z. The (111) surface has the lowest while the (110) exhibits 
the highest surface energy, as anticipated. On the other hand, the bulk 
energy e^ (the lowest curve) displays an opposite trend (i.e., e q 
increases with Z). This curve reflects structural characteristics of an 
fee crystal and depends only on the atomic configurations of the system 
.15]. 
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Calculated results for the multilayer relaxation of (110), (100) and 
(111) index planes are shown in Figs. 2a, 2b and 2c, respectively, where 
percentage variations in the interlayer spacings d^> ^23 anc * ^34 were 
plotted versus Z values. The most dramatic relaxation took place for the 
(110) surface. For Z < 0.1, all three interlayer spacings of the (110) 
surface were found to be larger than the bulk value. On the other hand, 
for Z > 0.1, d^ and d^ become progressively smaller than the bulk 
interlayer spacing, while the d£^ continues to expand (see Fig. 2a). The 
d^ exhibits the largest variation among all. A relatively smaller relax- 
ation took place for the (100) surface. Up to Z z 0.2 all three interlayer 
spacings remain larger than the corresponding bulk value. For Z > 0.2, 
the d^ exhibits some contraction, while the variations in the d^^ and 
d^ with respect to the bulk value remain negligibly small (see Fig. 2b). 

For the (111) surface, we found the smallest relaxation; however, the vari- 
ation in the d^ still displays a similar trend as in (110) and (100) 
cases. Up to Z ~ 0.3, again all three interlayer spacings were found to be 
somewhat larger than the bulk value. In this case, percentage variations 
in d 23 and d^ for z > 0.3 are negligible. Only the d^ exhibits a 
small contraction for Z values larger than 0.3 (see Fig. 2c). 

The quantity 9 which is defined as | a/e | is a useful property and 
can be regarded as the "relative" surface energy [36]. Figure 3 shows the 
values of 9 plotted versus Z for (110), (100) and (111) surface planes. 
The overall trend is similar to the variation of a in Fig. 1. The (111) 
plane remains as the energetically most favorable surface, consistently; 
and (110) is the least stable one bearing the highest surface energy. 

Values of 9 calculated for Z = 0 (which is a special case in our 
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approach) are in good agreement with results reported by Benson and Claxton 
[36]. 

Results obtained in this multilayer relaxation calculation (for the 
low index planes of fee crystals) are in good qualitative agreement with 
recent experimental findings [23-31]. For the Z = 0.3 case, the calcu- 
lated percentage variations in dj^» d^^ and ^34 f° r t * ie (HO) plane 
were found to be -9.98, +5.52 and -1.81, respectively; and, for the (100) 
plane, -1.63, -0.00 and +0.02, respectively. For comparison, Davis and 
Noonan [23] using a LEED technique obtained Ad^ = “10.0 + 2.5% and 
^23 = * 2*5% for the (110) surface; and Ad^ = “1*1 * 0.4% and 

^23 = * 0.6% for the (100) plane of copper. In another study based 

on a High Energy Ion Scattering experiment, I. Stensgaard et al. [27] 
obtained Ad^ “ -5.3 * 1.6% and ^ 22 , = +3.3 ± 1.6% for the Cu (110) 
surface. Furthermore, studies of the Cu (110) surface, which have been 
carried out by Adams et al. [30] based on LEED measurements, produced 
Adj^ = -8.5% and ^23 = ^»3%. general, we found a very small surface 
relaxation for the (111) plane. This conclusion is well supported by 
various experimental reports on the (111) surfaces [29,27]. Another 
comparison with experiments is possible for the surface energy of copper. 
The experimental value of 0 (surface energy per atom/absolute cohesive 
energy per atom) is 0.172 [37]. Our multilayer relaxation calculations for 
Z = 0.3 produced 0 = 0.185 which can be considered in fair agreement 
with the above results for the ( 111 ) surface plane. 

These comparisons indicate that, for the present potential energy 
function, the best Z-value for Cu should be around 0.3. Calculated 
- results also agree with the more recent experimental work by Andersen et 
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al. [29J for the multilayer relaxation of AA (110) surface. Furthermore, 
our calculations support the multilayer relaxation results obtained by 
Landman et al. using an electronic charge density method based on more 
accurate quantum mechanical considerations. There, they also reached the 
same conclusion: that the central two-body potentials alone are not suffi- 

cient to describe the surface region properly [33,34], 

However, care must be exercised in comparing these calculated multi- 
layer relaxation results with experimental findings. The calculations of 
the d„ 's and of 9 require only the knowledge of m , n and Z param- 
eters. If the absolute values of e,, e„ and d.. are desired, e and 

A O l j 

Tq are also needed. All these potential energy parameters must be deter- 
mined from the physical properties of the materials under consideration. 

The relationship between these potential energy parameters and the 
crystal structure has been shown elsewhere [5], The stability region of a 
crystal structure depends on the values of m and n as well as Z. In 
the present investigation we analyzed the dependence of the multilayer 
relaxation on the value of Z, the three-body intensity parameter. How- 
ever, we believe that an analysis of the effect of m and n values (of 
the two-body potential) would be an interesting study. 

Conclusions 

The most significant outcome of the present investigation is the 
understanding of the important effect exerted by the many-body forces on 
the multilayer relaxations of surfaces. Consideration of the three-body 
interactions corrects the shortcomings of pair potential models, and furn- 
ishes relaxed surface configurations with varying interlayer spacings which 
agree (at least qualitatively) with various experimental measurements. 
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For the low index surface planes of an fee crystal, the importance of 
the multilayer relaxation increases as (111) -*■ (100) -» (110); i.e., with 
decreasing atomic density of the plane. The first interlayer spacing d^ 
is found to be increasing for smaller three-body interaction intensities 
(Z), but decreasing for systems with higher Z. The second layer spacing 
622 tends to increase monotonically with increasing Z on the (110) 
surface, but decreases slightly on the other two surfaces. The percentage 
variation in the d^ is generally negligible, except for the (110) case 
where it is noticeable and decreases with Z. 

The potential energy employed in this study has a relatively simple 
functional form; therefore, it may be used in lengthy computer simulation 
calculations related to surfaces. 
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Figure 


Fig. 1 


Fig. 2 


Fig. 3 


Captions 

Change in energy per atom versus Z for different surfaces and at 
bulk. Both quantities, energy and Z, are in reduced units. 

Percentage change in interlayer spacings for different surfaces as 
a function of Z. Thick solid line, thin solid line, and dashed 
line represent percentage changes in d^* ^23* an< * ^34 
respectively. 

Ratio of the energy per atom at the surface to the absolute bulk 
value, 9, for different surface planes as a function of Z. 


24 



ngy 






CHAPTER IV 


Evaluation of Parameters 

In order to calculate the total potential energy for a system of N 
particles in a given configuration, the parameters (e, r^, m, n and Z) of 
eqs. ( 2)— ( 5) must be known. By definition, these parameters are materials 
constants and depend on the atomic species involved in the interaction. 

For a crystalline system the total potential energy expressed by eq. (6) 
may be simplified employing Lattice sum representations: 


N_ E 
2! m-n 



r r 0' ) n 


a 1 + 2® T 

n J d 9 \ 


( 11 ) 


where d is the nearest neighbor distance in the crystal. The lattice sum 
A's and are given by 
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r. . - r., 
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These lattice sums are simply geometrical factors and only depend on 
the crystal structure [5,7], Table 1 tabulates calculated values of- 
for different JL, along with the values for various crystal 
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geometries. In the evaluation of parameters we also used the stability 
criterion which is given by: 


'■av J T 


( 14 ) 


where V denotes the total volume of the system. This relation is exact 
only at T = 0°K. However, it has been always assumed that at relatively 
lower temperatures the minimum of the $ with respect to V coincides 
with the minimum of the free energy. Furthermore, the second derivative of 
the energy is related to the bulk modulus by: 

B - V(^f) . (15) 

6V 

For simplicity as well as for proper comparison of the results with 
other Leonard- Jones ion calculations, in this project, generally, we con- 
sidered that the exponents in the two-body potential function are m = 12 
and n = 6. Now, with this assumption, the remaining three unknown param- 
eters can be calculated from eqs. (11) through (15), using experimental 
cohesive energy and the bulk modulus. However, in some cases, these non- 
linear simultaneous equations cannot be solved. Therefore, for those 
cases, we used either surface energy values or experimental small cluster 
data for the fitting procedure. 
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Table 1 


Lattice Sums of the Mie Potential and the Axilrod-Teller 
Potential for Various Structures (One Component System) 



HCP 

FCC 

BCC 

SC 

DIA 

GRAH 

A 4 

23.616 

25.5946 

21.1685 

15.485 

9.5795 

5.4351 


16.883 

16.8807 

14.6913 

10.333 

6.2862 

3.8825 


14.449 

14.4481 

12.2495 

8.3994 

5.1153 

3.3895 

*7 

13.360 

13.3590 

11.0539 

7.4669 

4.5944 

3.1910 

*8 

12.803 

12.8019 

10.3551 

6.9458 

4.3310 

3.0993 

*9 

12.493 

12.4926 

9.8945 

6.6289 

4.1904 

3.0534 

A 10 

12.312 

12.3113 

9.5644 

6.4261 

4.1110 

3.0294 

A 11 

12.201 

12.2009 

9.3132 

6.2923 

4.0655 

3.0164 

A 12 

12.132 

12.1319 

9.1141 

6.2021 

4.0389 

3.0092 

A 13 

12.088 

12.0877 

8.9518 

6.1406 

4.0233 

3.0052 

A 14 

12.059 

12.0590 

8.8167 

6.0982 

4.0140 

3.0030 

A 15 

12.040 

12.0400 

8.7030 

6.0688 

4.0086 

3.0017 

A 16 

12.027 

12.0274 

8.6063 

6.0483 

4.0051 

3.0010 

A 17 

12.019 

12.0198 

8.5236 

6.0339 

4.0031 

3.0005 

A 18 

12.013 

12.0130 

8.4525 

6.0239 

4.0019 

3.0003 

A 19 

12.009 

12.0094 

8.3914 

6.0168 

4.0011 

3.0002 

^0 

12.006 

12.0063 

8.3386 

6.0119 

4.0007 

3.0001 

a 

19.175 

19.1697 

14.7719 

6.6138 

1.6647 

0.1010 


7200 

6912 

8192 

8000 

5832 

5040 


Indicates the number of atoms which were considered to compute 
the lattice sums. 


For covalent materials such as Si, C and SiC we calculated the 
parameters using eqs. (11) and (14), plus the small cluster data. The 
calculated results are given in Table 2. In the rest of this report, those 
parameters will be employed in various simulation calculations. 
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Table 2 


Constant Set in the Two- and Three-body Parts 
for Si, C and SiC 


m(Si,Si) 

= 12 

n(Si,Si) = 6 

£ Si,Si = 3.228(eV) 

RSi, Si = 2.29505(A) 

m(Si,C) 

- 12 

n(Si,C) = 6 

e Si,C = 4 ‘ 462(eV) 

Tq 1 ’ C = 1.74(A) 

m(C,C) 

- 12 

n(C,C) = 6 

c QC = 6.2295(eV) 

Tq’ 0 = 1.4806(A) 


2 Sl,Si,Si - 3991. 7 (e v • A 9 ) Z slsljC - 800.0(eV • A 9 ) 

Z Si c c = 300.0(eV . A 9 ) Z Q CQ = 191.73(eV • A 9 ) 


To investigate the applicability of these calculated parameters, 
several tests were carried out. 

Under normal conditions, silicon has a diamond cubic structure. Under 
high pressures, however, it undergoes to a phase transformation. The high 
pressure form of silicon has a tetragonal p-tin structure. This struc- 
tural transformation is accompanied by a large volume decrease (~22.7%). 
First, a simulation calculation was performed to analyze this phase trans- 
formation phenomena. For silicon we considered four different crystalline 
structures (including diamond cubic, p-tin, fee and bcc) and using the same 
set of parameters we calculated the total energy as a function of the total 
volume. This is shown in Fig. 1 which clearly indicates the energetically 
most favorable structure is the diamond cubic. Under higher pressures 
(i.e., for smaller volumes), the p-tin structure becomes energetically 
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more stable. According to this calculation, other structures may become 
stable at much higher pressures. The calculated change in volume for the 
transformation is in good agreement with the experimental data. 

Next, the melting process of silicon was simulated using again the 
same set of parameters. A constant pressure Monte-Carlo calculation was 
employed to simulate a silicon system which contained 64 atoms with 
periodic boundary conditions in all three directions. For every incre- 
mental temperature step (from 300°-5000°k), the system was equilibrated 
which was monitored by the variation of the total energy. In Fig. 2, the 
total energy, relative volume and the bulk modulus values are plotted as a 
function of the temperature. All three of these functions exhibit changes 
of slope around T = 2000°K. Furthermore, calculated radial distribution 
functions as well as trajectory plots of the particles indicate that the 
melting takes place at ~2000°K. This is somewhat higher than the experi- 
mental melting point of silicon, however, the ability for the potential 
energy function to provide a proper volume decrease during the melting 
process was considered as an important accomplishment. 

Further checks for the adequacy of the parameter set involved also the 
investigation with pure carbon. The potential energy function for C with 
parameters tabulated in Table 1 can provide two closely spaced energy 
levels for diamond and graphite. Calculated energies are somewhat larger 
than the experimental values; however, the calculation predicts that the 
lower lying state belongs to graphite. Also, for the case of SiC the 
potential energy provides two almost degenerate energy states corresponding 
to a and {3 forms of SiC, consistent with experimental data. 
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Figure Captions 


Fig. 


Fig. 


Total energy curves for the four structures of Si as a function 
of the atomic volume. Dashed line is the common tangent of the 
energy curves for the diamond cubic and p-tin structure (c/a = 
0 . 546 ). 


(a) The total energy, (b) the volume and (c) the bulk modulus of 
the system as a function of temperature at the equilibrium state. 
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CHAPTER V 


Simulation for Surface Properties 

Detailed atomistic level information about clean solid surfces is very 
important in the analysis of various surface processes. Calculation of 
surface properties from interatomic potentials is an involved procedure due 
to the fact that the role played by many-body forces at the surface region 
is not negligible. For an accurate simulation calculation these many-body 
effects must be accounted for properly [35]. In particular, for substances 
like Si or SiC which were considered in this investigation, the many- 
body interactions are shown to be very significant. 

Surfaces were generated in the computer as abrupt discontinuities of 
the crystalline bulk. To obtain relaxed (or reconstructed) surfaces, this 
initial configuration was permitted to relax under the Monte-Carlo code. 

In general, we applied periodic boundary conditions in two-directions 
which provided an effectively infinite exposed surface in the desired 
direction. Relaxation procedures were carried out for the top three to 
five atomic layers (depending on the surface geometry). The rest of the 
atoms in the system were fixed in their original lattice sites. However, 
the fixed atoms contributed fully to the total energy calculation. 

Si (100) Surface: 

Calculations for Si(100) surface were performed at two different 
temperatures employing a Monte-Carlo procedure. At T = 298°K the equi- 
librated structure of Si(100) exhibited a c(2 x 2) reconstruction. As 
it is shown in Fig. 1, the atoms at the top layer tend to form pairs. This 
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feature is quite similar to reconstructed models of Si(100) suggested by 
various investigators [38-41]. It is now accepted that the general 
tendency in the structural reconstruction of the Si(100) surface is the 
formation of dimers as originally suggested by Schlier and Farmworth [42]. 
More than one type of reconstruction patterns may coexist on Si(100) 
surface (e.g., not only (2 x 1) } but possibly c(4 x 2), p(2 x 2) and 
c(2 x 2) superstructures as well). Our result c(2 x 2) represents one 
bf these reconstruction models. 

Atoms located in the second and third layers remained more or less 
stationary. Interplanar spacings between the first and the second layers 
displayed a contraction, while other interlayer spacings were affected only 
marginally. Our calculated results for the interplanar relaxation are 
quite consistent with results reported by Yin and Cohen [43] who employed 
an ab_ initio self-consistent pseudopotential method. 

The surface energy calculated from an equilibrated Si(100) surface 

2 

was found to be 1386 ergs /cm which is in good agreement with experimen- 
tal data [37]. At low temperatures (T ~ 1°K) on the other hand, it was 
found that the surfce structure after relaxation preserved its (1 x 1) 
symmetry (i.e., no dimer formation was observed). However, the first 
interlayer spacing for the equilibrated Si(100) surface exhibited a 
contraction as in the previous case explained above. 

In the second part of this chapter, statics calculations were per- 
formed for the stress tensor components as well as the energetics and the 
structure for the (111) surfaces of Si and SiC. 

The major contribution to the stresses in a system of interacting 
particles is the strain derivative of the potential energy. Since the PEF 
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used is a sum of terms depending on only a few atoms at a time, the stress 
is also. It is thus possible to define a potential and stress for each 
atom in units of energy/atom. In the next section, these atomic quantities 
are plotted using the corresponding atomic positions to present a qualita- 
tive description of the highly non-homogeneous behavior of the surface or 
defect. 

A variety of excess energies and formation energies are reported in 
the following section. They have been calculated using the general expres- 
sion: 


E 


ex 

a 



(16) 


E^, is the total energy for the simulation of N^, atoms with bulk energy 
per atom <j>. C is a normalization factor, such as area or length. The 
final term in eq. (16) subtracts all other appropriate total excess 
energies. Thus, for example, the ledge excess, k, would be: 



(E t - N t <)> - A T y) 


(17) 


where y is the flat surface excess energy determined by a previous simu- 
lation, Aj is the total exposed surface area in the simulation with a 
ledge crossing it of total length L T . Of course, these excess energies 
are not free energies, but internal energies. All of the results of the 
next section are thus subject to modified interpretations depending on the 
possible entropic effects for the "real" system. 
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Silicon (111) Surface Properties: 

1) Perfect Surface 

When the diamond cubic lattice is terminated on a (111) plane, the 
unrelaxed geometry is as shown in Fig. 2(a) and 2(b). The pairs of layers 
(C-a), (A-b), or (B-c) are best considered as one "puckered" layer of 
"upper" and "lower” atoms, designated by open and filled circles respec- 
tively. The intra- and inter-layer spacings for the bulk are 6 = .7839 
and A = 2.3516 A. The region enclosed by dashed lines in Fig. 2(b) is 
the primitive two-dimensional surface unit cell. 

If we take these perfect crystal positions as the initial positions, 
then our relaxation calculations show a 27% contractin of the outermost 
layer to 6^ = 0.574, and a slight expansion of the other distances Aj , 

6^, by 1.6, 1.4, and 0.05 percent respectively. The layers remain flat 

and otherwise unreconstructed. Our calculation relaxed four puckered 

2 

layers of 50 atoms each (A^ = 319.25 A ) to a final total energy of 

2 

-1230.50 eV, yielding y = 1169 ergs/cm . Figure 2(c) shows the potential 
energy distribution in the first four surface layers; Although the surface 
energy is positive as it must be, notice that the potential energy of the 
"lower" atom s in the first layer is actually more negative (more strongly 
bonded) than atoms in the bulk due to the surface reconstruction. 

Figure 2(d) shows the stress, t^, for the first four surface layers, 
is the normal stress along the (Oil) direction (in the plane of the 
surface), and was found to be identical to the normal stress along the 
(211) direction. The surface is evidently in compression relative to the 
bulk. If the total stress in the first two layers is assumed to be dis- 
tributed over the region indicated in the figure, we can convert eV/atom to 


39 



energy/volume , with the result that T ]1 = ^1437 atm. We can draw two 
immediate conclusions: (i) if the crystal is expanded for any reason 

(thermal strain, substrate misfit, etc.), the surface excess energy will 
decrease, and (ii) defects which generate local tensile fields will be 
attracted to the surface. 

2) Perfect Ledges: 

There are two unequivalent high-density ledges which terminate 
terrace layers on the (111) surface: the (211) and (211). If we con- 

sider Fig. 1(b) as a terrace, then these ledges are to the left and right 
respectively. The dashed lines indicate the "broken bonds” along the 
ledge. We have simulated both of these ledges separately. The periodic 
boundary conditions were established such that the ledges were widely 
spaced at 26.6 A (8 surface unit cells). After relaxation, the ledge 

excess energies were found to be X _ = 0.30 eV/A and X — - 0.16 - 

( 211 ) ( 211 ) 

eV/A. 

One might initially expect that the (211) ledge would have a lower 
ledge excess energy due to the fewer number of "broken bonds." However, 
the (211) ledge undergoes a major reconstruction, in which alternative 
pairs of ledge atoms form dimers as shown in Fig. 3(a), with a spacing of 
2.33 A. It is this reconstruction which causes the (211) ledge to be 
energetically favored over the (211). This result agrees qualitatively 
with other theoretical calculations [44] and with cleavage experiments 
[45]. 
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The energy and stress distributions around the lower energy (2ll) 
ledge are shown in Fig. 3(c)-(e). The plots are of the quantities in 
eV/atom for all the atoms on one (111) layer; starting from the right of 
Fig. 3(b) at a ledge, moving left across the terrace, moving under the next 
ledge, etc. There is a one-to-one correspondence between the positions of 
the atoms in the indicated layer in Fig. 3(b) and the positions of the 
peaks in the plots below it. 

The potential energy distribution, Fig. 3(c), shows only a small 
transition region as the layer passes under the ledges. The values for the 
energy per atom halfway between the ledges agrees to within 0.3% with the 
values from Fig. 2(c) for the perfect surface. From this, we might con- 
clude that the ledges are not interacting at this spacing. However, 
consider Fig. 3(d), the normal stress directed parallel to the ledge, 
x . The dotted line indicates the value for the larger peak from Fig. 
2(d). The ledges are thus interacting even at this large spacing, since 
the stress does not achieve the perfect surface value. The ledges act to 
decrease the compressive character of the upper surface in the ledge- 
parallel direction; however, the overall effect is not as simply character- 
ized. The most distinctive effect due to the ledge is the development of 

the shear stress field, x , surrounding it, as seen in Fig. 3(e). 

xy 

Several features should be emphasized. First, the reconstruction of 
the ledge opens up large holes along it which may become favorable sites 
for impurity adsorption. Second, the stresses, and in particular the shear 
stress, indicate how the ledge may affect other defects in the system. 
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3) Kink Site: 


Investigation of the behavior of kink sites on the (211) ledge 
has just begun. The kink site excess energy is 1.67 eV and it has only a 
minor effect on the stress distributions around the ledge. It is important 
to note that the '’motion'' of the kink site requires the addition of four 
atoms, due to the dimer reconstruction of the ledge. The activation energy 
for Si adatom attachment at the kink site will be strongly influenced by 
this fact. 

4) Adatoms and Surface Vacancies: 

The two most likely adsorption sites for single Si adatoms on the 
otherwise flat (111) surface are the three-fold coordinated sites [46]. 

The "cradle" site has a "lower" atom directly underneath it while the 
"hole" site does not (see Fig. 4). The surface vacancy is created by 
removing one "upper" atom. The defect formation energy E^, as calculated 
by eq. (16), can be considered as the energy change in the system using a 
kink site (i.e., the bulk potential energy per atom) as a source or sink of 
atoms, and is a measure of the stability of the surface to the formation of 
these point defects. 

The formation energies were determined for an effective surface defect 
concentration of 2% (1 defect every 25 surface cells). The cradle, hole, 
and vacancy energies are, respectively, 0.84, -0.31, and -1.19 eV. The 
perfect surface thus lowers its energy by 1.50 eV in creating a hole/ 
vacancy pair, and no adjacent ledge is required as a source or sink of, 
atoms. A further calculation showed a reduction to = -1.82 eV if the 
hole and vacancy associate on adjacent surface sites (surface Frenkel 
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pair). Our potential thus predicts that the (111) surface will spontane- 
ously roughen, or reconstruct, to some equilibrium geometry of adatoms and 
vacancies. This behavior provides a basis for the interpretation of the 
(7 x 7) reconstruction pattern observed on the Si (111) surface by LEED 
and STM experiments [46-48], Further analysis of this effect is in 
progress. 

The surprising stability of the surface vacancy depends on two 
effects. First, the vacancy generates tensile fields in the surrounding 
lattice and is thus favored in the highly compressive region at the 
surface. Even for the small concentration of defects used above, there was 
an 85% reduction in the surface compression (relative to the flat surface) 
for the vacancy, as compared to only a 6% increase in the surface compres- 
sion for a hold adatom. The second effect is a strong reconstruction 
around the vacancy as shown in Fig. 3. The lower atoms surrounding the 
vacant site pull together in a configuration similar to the dimer recon- 
struction on the (211) ledge. If we prohibit this reconstruction, the 
vacancy formation energy becomes +2.68 eV, which is consistent with other 
calculations for the unreconstructed surface vacancy energy [49]. 

B) Silicon Carbide (111) Surface Properties 

Many of the calculations described above have also been performed for 
the equivalent geometries on the (111) faces of SiC in the zincblend 
structure. We will briefly summarize the results obtained. 

The zincblend binary structure has two unequivalent (111) faces, which 
in this case can be referred to as C-rich and Si-rich. The latter under- 

2 

goes only a 5% contraction of the first layer, leading to y = 2544 ergs/cm 
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and t = 20170 atm. The C-rich face, however, develops an outer layer 

2 - 

contraction of 35% leading to y = 344 ergs/cm and x = -100619 atm. 

The crystal growth properties of these two surfaces should thus be quite 
different. This places enormous importance upon the initial nucleation 
even for the SiC since a given heterogeneous substrate may favor one 
orientation over the other. 

Since the C-rich face is so much more stable, we have concentrated on 
it. The (211) and (211). ledges on this face have excess energies of 
0.50 and 0.96 eV/A respectively. Analysis of the geometry of the 
(211) ledge shows that the dimer reconstruction is very minor compared 
to the silicon case, thus the (211) ledge is more stable as expected from 
the simple broken-bond argument. The stress distributions around the 
(211) ledge are similar to the silicon (2ll) ledge as in Fig. 3(c-e), 
except that the sign of the shear stress is reversed at all points. 

The formation energy for carbon and silicon adatoms in the hole sites 
on the C-rich face are -6.11 and +2.66 eV respectively, based upon the 
individual bulk energies per atom in the SiC crystal of <)> = -9.995 eV 

and 4> c = -7.6627 eV. The individual formation energies are somewhat 
arbitrary, but they do indicate that the carbon is much more strongly 
adsorbed onto the surface that the Si. The "smaller" carbon adatom is 
drawn down almost to the same level as the other C atoms on the surface, 
while the Si adatom remains approximately 1.3A above the surface. 

Very useful qualitative insights concerning surface energetics and 
surface processes can be gained from computer simulation studies using 
two-body plus three-body PEF's. The behavior of the stress distributions 
•around these same structural elements has been described. This is an 


44 



especially useful tool for describing and predicting the overall growth 
process in terms of the cooperative interactions of these basic structural 
elements. 

Several other features should be highlighted: 1) The Si (111) 

surface is in compression. 2) The SiC (111) C-rich face is in compression 
while the Si-rich face is in tension. 3) The Si (111) surface should 
reconstruct to some equilibrium configuration of surface Frenkel pairs. 

4) The (211) ledge on the SiC (111) C-rich surface is the most stable. 

5) The (211) ledge on the Si (111) surface is the most stable due to 
dimer reconstruction. 6) Both ledges develop well-defined shear stress 
fields underneath them. 
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Figure Captions 


Fig. 1 Relaxed configuration of the Si (100) surface calculated at 
T = 298°K. 

Fig. 2 (a) Side, and (b) top view of (111) surface of silicon. 

(c) Potential, and (d) surface stress distributions for layers as 
labeled in (a). 

Fig. 3 Properties of (2ll) ledge on Si (111) surface, (a) Dimer 

reconstruction, (b) Layer of atoms whose (c) potential, and (d-e) 
stresses are shown directly below. All quantities are in units of 
(eV/ atom) . 

Fig. 4 Relaxed equilibrium geometry of surface vacancy. Also shown are 
the cradle (^^) and hole (^ ) adatom adsorption sites. 
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Top view of the unrelaxed Si(100) surface 
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Top view of the relaxed Si(100) surface 
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FIGURE 4 
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CHAPTER VT 


Atomistic Modeling for Slip Formation 

The mechanical behavior of materials can be specified by macroscopic 
theories on the basis of a few material constants that provide an accurate 
description for the responses of materials to forces. However, these 
theories, for example, do not provide a microscopic level of understanding 
of the basic mechanisms involved in plastic deformation. In general, the 
theories of strength of materials, elasticity and plasticity lose much of 
their power when the structure of a material becomes an important consider- 
ation and the material can no longer be considered a homogeneous medium 
[50]. 

The relationship between mechanical behavior and microscopic structure 
of materials is very important. When mechanical behavior is understood in 
terms of microscopic processes, it is often possible to improve the mechan- 
ical properties of a material. A microscopic description of materials, as 
opposed to macroscopic theories, cannot be adequately defined by using a 
few material constants. Instead, the system must be described in terms of 
interatomic forces and the coordinates of the particles which constitute 
the material. 

Several studies based on the computer simulation of the mechanical 
behavior of materials exist. These studies have provided a better under- 
standing of various mechanisms (diffusion, crack propagation, dislocation 
motion and plastic flow) [51,52] at the atomistic level. 

Material failure is generally caused by fracture which follows yield- 
ing or plastic deformation. Slip is the simplest and the most common 
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example of plastic flow encountered prior to the ductile fracture of 
materials. In this investigation the effect of a uniaxial load exerted on 
a two-dimensional microscopic slab was analyzed. 

The system contained 400 or more particles. Each, particle was treated 
discretely. All particle neighbors up to 3.5 Tq were included in the 
energy force calculations. This procedure ensures that all neighbors up to 
the f if th-nearest neighbor will be included in the calculations. In the 
calculation procedure, periodic boundary conditions were not employed. 

This allowed us to examine the surface region reconstruction and the 
formation of edge dislocations during the elongation process. 

Initially the system was generated in a rectangular shape in its 
equilibrium configuration. Then the system was elongated in a step-wise 
fashion by imposing a uniaxial load in small increments. After each incre- 
mental elongation, the system was allowed to equilibrate fully during a 
relaxation period. All the mobile particles were fully equilibrated after 
every incremental elongation by using a force minimization technique. The 
force acting on each particle was calculated; the particle was then moved 
in the force direction until the resultant force acting on it became 
virtually zero. This procedure was repeated sequentially for each moving 
particle up to the complete equilibration. 

The discrete model used in this investigation produced results that 
are consistent with those from macroscopic theory. The results herein 
indicate that slips, which are known to be the simplest type of plastic 
deformation in crystalline bodies, occur predominantly on rows with a 
higher density of atoms (along the close-packed rows). In addition, the 
calculations for perfect two-dimensional triangular crystals with uniaxial 
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loads imposed in the [Oil] , [112] and [514] directions showed that 
the rule of maximum resolved shear was observed [50]. These simulation 
results as expected are a confirmation of the macroscopic theories and 
illustrate the involvement of dislocations in the slip formation process. 

In all cases studied, the system with point defects experienced slip 
formation at smaller strains than the corresponding perfect crystals. 
Vacancies located near the surface moved to the surface before the slip 
occurred. However, vacancies in the interior regions moved to the surface 
during the slip process. 

The details of this static simulation calculation were published 
recently (see reference [53]). 

To analyze the process of slip formation in "real time," in addition 
to the static approach described above, we also made use of a molecular 
dynamics technique. A similar model was taken into consideration and the 
tensile load on the system was generated in small incremental elongations. 
However, in this case, the system was equilibrated with a dynamic code at a 
finite temperature. The results obtained in this investigation were found 
to be basically similar to those obtained by the static method (except for 
a small temperature effect). This dynamic simulation provided some addit- 
ional information about the involvement of dislocations in the slip forma- 
tion process and its time dependent characteristics. 
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CHAPTER VII 


Molecular Dynamics Calculations for Crack Propagation 

In this investigation we analyzed the process of crack propagation in 
two-dimensional lattices using a molecular dynamics technique. Simulation 
calculations were carried out considering systems containing approximately 
2400 discrete particles interacting via two-body potentials. In addition 
to energetics, forces and instantaneous position vectors, we also calcu- 
lated stress components for every particle in the system. General char- 
acteristics of the model considered in this study are basically similar to 
the one used in the previous chapter for the slip formation calculations. 

Simulations for two-dimensional systems are relatively easier to 
analyze than results for three-dimensional systems. First of all, 2D 
systems contain a smaller number of particles and, therefore, require less 
computer time. Results can be represented by simple 2D plots and problems 
arising due to the multi-particle character of the system are easily 
identifiable. Thus, the analysis of the 2D system provides considerable 
information not only about the microscopic nature of the crack growth 
phenomenon, but also provides some knowledge about "how to interpret the 
simulation results." The question of the credibility of these 2D results, 
of course, remains unanswered. At this stage, it is not known how to 
extrapolate results obtained from a 2D simulation to a 3D domain. However, 
the results obtained in this study together with several other reports 
[54,55] in the literature indicate that 2D systems, in most cases, do 
exhibit characteristics similar to 3D systems. 
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Atomistic level analysis of the crack propagation process using 
computer simulation techniques has been the subject of several earlier 
investigations. In the literature we could find only a few reports 
relevant to the study carried out in this investigation. In the report by 
Ashurst and Hoover [54] , the fracture phenomenon was investigated based on 
a truncated Hook's law force. They have found that, even with this very 
simplistic force law, their static simulation furnished results for energy, 
entropy, stress concentration and crack structure all to be consistent with 
expectations from macroscopic elastic theory. 

The other relevant and more recent study was reported by Dienes and 
Paskin [55]. In this modeling study they also considered a 2D triangular 
lattice with particles interacting via the Lennard-Jones function. A crack 
has been introduced in the interior of a pre-stressed sample. The crack 
was initiated by "cutting" the bonds between a given number of atoms at the 
central portion of the sample. The interatomic potential was artifically 
set to zero between these atoms. According to their report, the condition 
would correspond to the insertion of a. very thin knife to create the 
crack. Furthermore, in the energy and force calculations, they only 
considered nearest neighbor interactions (by taking R cut “1*6 r o^* 
their model, the crack was aligned parallel to close-packed rows and 
displayed a linear path in its propagation. Finally, they found that their 
results are quantitatively good at the early stages of the propagation 
process. 

The main objective of this study is to investigate the crack propaga- 
tion phenomenon at an atomistic level to understand and resolve various 
mechanism involved in a crack tip process. In the first part of this 
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simulation study (to test the system under consideration) we calculated the 
stress-strain characteristics of a perfect two-dimensional lattice at a 
finite temperature. This calculation was carried out primarily for 
comparison with the results obtained in the second part of the study where 
calculations were performed for systems with pre-existing cracks. 

(i) Perfect Lattice: 

As a perfect 2D lattice, the basal plane of an hep lattice was 
taken into consideration. A system of 2400 particles in a rectangular 
shape (80 x 30) was first generated in static equlibrium. A tensile load 
was applied in the [112] direction, which is the close-packed direction. 
This direction was also chosen as the x-direction in our cartesian 
coordinate system. The load was imposed in small incremental strains (in 
this case elongations) of 0.01. This was performed uniformly throughout 
the system by factorizing all the x-components of the position vectors 
describing the system. In the x-direction, periodic boundary conditions 
(PBC) were applied to provide continuity for the system (in the tensile 
direction), and also to furnish two free surfaces in the y-direction. In 
a general sense, the imposed PBC provides the desired tensile strain on the 
system. 

The system was relaxed after every incremental strain by a molecular 

dynamics code. A cut-off radius, R _ , of 2.86 r n was considered for 

cut 0 

the energy and force calculations. This R ^ is between the fourth and 

cut 

fifth neighboring shells surrounding the central atom and provides approxi- 
mately 30 neighbors. The reduced time step was taken as 0.01 and the 
reduced temperature was T* = 0.02 (to compare with real systems; e.g. , for 
copper these represent 2.5E-15 seconds and 100°K, respectively). 
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For every strain value, the stress components of particles were 
estimated as derivatives of the potential energy per particle. The 
interaction energy for a particle i was calculated as: 

M. 

4> i = [ u ( r ij) < 18 ) 

where M is the total number of neighbors of atom i and r„ denotes 
the distance between particles i and j. To represent the pair inter- 
actions we employed the Lennard-Jones function: 

u(r ) = t[(^-) 12 - 2(^-) 6 ] (19) 

J ij r ij 

with e and r^ denoting the energy and the internuclear distance, 
respectively, at equilibrium. 

For each particle the stress components were calculated considering 
Lagrange strain parameters. For example, the stress component, for a 
particle i, in the x-direction is given by: 


1 

a 

xx 


, M r n 

- f I K^ 2 -) 

a n -i j .2 


12 


ij ij 


( 20 ) 


where a^ denotes the area per particle and x^ is the x-component of 
the position vector for the particle i. 

The calculated stress-strain curve for the perfect lattice case is 
shown in Fig. 1, up to e = 0.09. 
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(ii) 2D Lattice with an Existing Crack 

A lattice with an initial surface crack was generated by removing 
9 particles from the surface region of a perfect lattice (see Fig. 2). 

This system, now with an existing surface crack, was elongated and relaxed 
by the molecular dynamics code in a similar way explained above for the 
perfect case. First, the effect of the temperature on the stress-strain 
curve was analyzed up to e = 0.03. Figure 3 shows two curves, dotted and 
solid, representing the stress-strain curves for T* = 0.1 and T* = 0.02, 
respectively. The shift in the dotted curve (high temperature curve) is 
mainly due to the thermal expansion. For lower strain values, these curves 
represent fully equilibrated systems. However, for strains higher than 
0.02, systems may require additional relaxation times to equilibrate 
completely. The difficulty involved in attaining the equilibrium is mainly 
due to large fluctuations displayed by the stress values calculated as 
derivatives of the energy. At this stage, we believe that the general 
trend exhibited by these two curves is sufficiently accurate for the 
present investigation. Any further incremental elongations (in addition to 
e = 0.03) cause the crack to propagate. Determination of the critical 
strain, i.e. , the strain at which the crack first starts to propagate, is 
difficult to assess. For this purpose we performed three separate runs 
with three different pre-s trained systems, namely with e = 0.03, e = 0.035 
and e = 0.04, all at T* = 0.02. The 2D lattice with the surface crack 
was strained in one single step from its original length up to 3.0, 3.5 and 
4.0% elongations. In the case of e = 0.03, the crack did not exhibit any 
growth and the overall configuration of the systems remained unchanged up 
to 3500 time steps. However, for both e = 0.035 and 0.04 cases, the 
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propagation of the crack took place. In these prestrained cases, we 
simulated the system under nonequilibrium isothermal conditions. For the 
e = 0.035 case, the crack growth first initiated after 1000 iteration 
steps. Figure 4 displays the stages of this relaxation process up to 3200 
iterations, at which the system reached almost to an equilibrium state. 

The darker circles in the figures represent particles with higher 
stresses. For the e = 0.04 case, on the other hand, the crack propagated 
much earlier (obviusly because of the high strain imposed initially). The 
crack started growing first at the 500th iteration step and the system 
attained an equilibrium state at approximately 2400 iteration. The stages 
of this propagation process are shown in Fig. 5. Again, the darker circles 
display particles with higher stresses. In both cases, the particles at 
the crack tip exhibited high stresses consistently. Furthermore, the crack 
propagated along the close-packed rows of the lattice and, at the same 
time, tried to remain perpendicular to the applied load by choosing a zig 
zag path. These behaviors are very much consistent with experiments and 
theories based on macroscopic considerations and, therefore, indicate the 
adequacy of the present atomistic level simulation procedure. The relaxa- 
tion of the system can be' followed in Fig. 6 where the average stress is 
plotted versus the iteration steps. The oscillatory behavior of this curve 
is a temperature effect mainly due to vibrational motions displayed by 
individual particles in the system. From the snapshops shown in Fig. 5, we 
also calculated the velocity of the crack propagation. The curve in Fig. 7 
represents the variation in the crack propagation velocity as a function of 
the calculated average stress. The upper range of this curve is near the 
velocity of sound propagation. This is expected according to a report by 
As hurst and Hoover [54]. 
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Figure Captions 


Fig. 1 Stress-strain curve for the perfect lattice. 

Fig. 2 Two dimensional lattice with an existing surface crack. 

Fig. 3 Temperature effect on the stress-strain curve. 

Fig. 4 Snapshots for the crack propagation with e = 0.035. 

Fig. 5 Snapshots for the crack propagation with e = 0.040. 

Fig. 6 Variation of the averaged stress as a function of time steps. 
Fig. 7 Crack propagation velocity plotted versus the averaged stress. 
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